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1. Introduction 

Variational and diffusion quantum Monte Carlo (VMC and DMC) are stochastic 
quantum Monte Carlo (QMC) methods for evaluating expectation values with many- 
body wave functions pp. The accuracy of VMC is quite limited in practice, but its more 
accurate and sophisticated cousin the DMC method is capable of giving very accurate 
results, often retrieving over 90% of the correlation energy. The computational costs 
of fermion VMC and DMC calculations scale approximately as the third power of the 
number of particles, making it feasible to deal with hundreds or even thousands of 
particles and allowing applications to electrons in condensed matter. 

DMC has been applied to a wide variety of condensed matter systems, including 
electron gases [21 El S] , nanocrystals [51 Ej, molecules on solid surfaces Ej, defects in 
semiconductors [9j [El E] , solid state structural phase transitions [12] , and equations 
of state [131 HS1 US]. DMC calculations have provided accurate benchmark results 
for many systems. QMC algorithms are intrinsically parallel and are ideal candidates 
for utilising the petascale computers which are now becoming available. Our CASINO 
code [17] has already achieved efficient parallelisation on machines with thousands of 
processors. 

Calculating energy derivatives such as atomic forces is very important in quantum 
mechanical simulations. Forces are used in relaxing structures, calculating their 
vibrational properties, and performing molecular dynamics simulations. It has, however, 
proved difficult to develop accurate and efficient methods for calculating atomic forces 
within DMC. Difficulties have arisen in obtaining accurate expressions which can readily 
be evaluated, and the statistical properties of the force expressions are less advantageous 
than those for the energy. In this paper we study approaches for calculating derivatives 
directly within QMC. Finite-difference energy calculations within QMC using correlated 
sampling methods [18] provide a different route to the same goal, and such methods have 
been developed by other workers [191 120] . 

2. The VMC and DMC methods 

In VMC the energy is calculated as the expectation value of the Hamiltonian with an 
approximate many-body trial wave function ^t- For a real the energy is 



where H is the many-body Hamiltonian and R denotes the vector of particle coordinates. 
To facilitate stochastic evaluation, Ey is written as 



Ey 



J^ T (R)H(R)^ T (R) dK 
/# T (R)# T (R) dR 



(1) 



Ey 



Jfr T (R) 2 E L (R) dR 
J^ T (R) 2 dR 



(2) 



where the the local energy is 



E L = V^HVt ■ 



(3) 
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The calculation proceeds by the generation of N configurations Rj sampled from the 
probability distribution 

^t(R)^t(R) m 
vl J /# T (R)tf T (R)dR' { ) 

using, for example, the Metropolis algorithm, and the energy is evaluated as 
1 N 

Ey^-^E^R,). (5) 

JV i=i 

Equation (j2J) with the Py of equation (jlj) is an importance sampling transformation 
of equation (prj which exhibits the zero variance property. As approaches an exact 
eigenf unction, Ei, becomes a smoother function of R and the number of configurations, 
N, required to achieve an accurate estimate of Ey is reduced. If \&t is exact, equation (j2J) 
gives the exact result even for a single configuration. Although this ideal limit cannot 
be reached in non-trivial calculations, it is expected that using a zero- variance estimator 
will lead to improved statistics. 

DMC reduces the systematic bias inherent in VMC by an evolution of the wave 
function in imaginary time. Such projector methods suffer from the infamous "fermion 
sign problem" in which the wave function decays rapidly towards the lower energy 
bosonic ground state. Stable fermionic behaviour may, however, be achieved by using 
the "fixed-node approximation" [2T) [22] in which the nodal surface of the DMC wave 
function $ is constrained to equal that of \&t- An importance sampling transformation 
is also used and the algorithm generates configurations distributed according to the 
"mixed" probability distribution, 

_ *t(R)$(R) , . 

Pd(R) - /* T (R)$(R)dR • (6) 

The DMC energy is given by 

_ J E L dR 

Eb ~ /<M>dR ' (?) 
which is evaluated as 

1 N 

En^-Y^E^)- (8) 

i=i 

Both the VMC and DMC energies are upper bounds on the exact ground state energy. 
The DMC energy is, however, more accurate as it is bounded from above by the VMC 
energy. 

We also calculate expectation values with the pure probability distribution, 

p (m _ <&(R)<&(R) (q) 

Pp(R) " ;$(R)$(R)rfR • (9) 

Pure expectation values can be obtained using several methods: the approximate (but 
often very accurate) extrapolation technique [IB], the future- walking technique [23] 
which is formally exact but statistically badly behaved for large systems and poor trial 
wave functions, and the reptation QMC technique [23], which is formally exact and 
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well behaved, but quite expensive. The extrapolation technique can be used for any 
operator, but the future-walking and reptation techniques are limited to spatially local 
operators. 

3. Forces in VMC 

The derivative of the VMC energy of equation ([2]) with respect to a parameter A in the 
Hamiltonian can be written as 



d3 



v 



1 dH 7 (H-E L )d^ T 

■Wt 1 



* T dA \Pt dA 



dR 



dA 



J ^ T ^ T dR 

Wt dA 



dR 



(10) 



This expression may be evaluated as an average over the probability distribution Py of 
equation (J4]). The term involving dH/dX is the Hellmann-Feynman force and the others 
are Pulay terms which contain the derivative of \&t with respect to the parameter in 
the Hamiltonian. Calculating VMC forces including the Pulay terms is straightforward, 
although it requires the evaluation of d\l/T/dA. 

Equation (flOl) can be written in a number of different ways which are equivalent 
when the expectation values are evaluated exactly. The particular form chosen satisfies 
the zero variance condition if the configurations are sampled from Py and H and 
dH/d\ are taken to act to the right. If and d^x/dA were exact, the exact 
force would be obtained from averaging over any number of configurations. It is 
expected that this property will reduce the statistical noise in practical calculations. 
Equation ( TlOl) has been used by a number of groups to evaluate forces within VMC 

[SSlESlETlEHlEHlEniET]. 

4. Forces in DMC 

The derivative of the DMC energy of equation (J7|) can be written as 



This expression includes the derivative of the DMC wave function d$ / dA which cannot 
be evaluated within the standard approach. A number of studies have used equation ffTTj) 
with the approximation [32] 

Id* 1 d^ T , . 
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which leads to an expression which is no more difficult to evaluate than the VMC one 
of equation ffTOl) . This approach gives reasonable results if \1/t and d^x/dA are accurate 
enough, but it leads to errors of first order in (\1/t — 3>) and (d^x/dA — d$/dA), which 
are often significant in practice. 

An alternative approach is to evaluate the forces using the pure distribution 
Pp of equation (Q, which is proportional to $$. It is more expensive to generate 
configurations distributed according to the pure distribution than the mixed one, but 
it brings significant advantages because less severe approximations are required in 
evaluating the forces. The energies evaluated with the mixed and pure distributions 
are equal, i.e., 

D Jty T <S>dR JQQdR ' 1 ' 

so that the DMC force evaluated with the pure distribution is 

dH „ t ,« „ , d$ 

(14) 



d^o /*dA*^, J^-^dA^ 



dA J dR J <£>$dR 

It turns out that the Pulay term in this expression can be written as an integral 
over the nodal surface of \l/x [131 [MJ [35] . This fact was first pointed out by Huang et al 
[33], and Schautz and Flad [34] provided an exact expression for the Pulay term, 

d$ 



1J 


f 1 
s 




d$ 

dX rfS 


2 


f $$ciR 



(15) 



where S denotes the nodal surface and V$ is the gradient of $ obtained on approaching 
the nodal surface from inside the nodal pocket along the direction normal to the surface. 
Unfortunately this expression cannot readily be evaluated because it involves (i) an 
integral over the nodal surface and (ii) the quantity d$/dA. 

One of the achievements of Badinski et al [35] was to show that d$/dA can be 
eliminated from equation ( Tl5|) because the nodal surfaces of $ and \l/x must be equal 
for all values of the parameter in the Hamiltonian. They obtained the result 

d$ f ^ |V* T | d^ T 7 c 



$$dR 2 / $$rfR 



which is the average over the pure distribution of a quantity written entirely in terms of 
^t- The difficulty of evaluating an integral over the nodal surface remains and Badinski 
et al [35| [31] attempted to circumvent this problem by writing the nodal surface integral 
as a volume integral. They noted that the second term on the right hand side of equation 
(fi~l| can also be written as a nodal surface integral, and that it is related to the pure 
distribution nodal term of equations (Tl5l) and (Tl~6l) by an extrapolation approximation. 
The standard extrapolation approximation is 

= 2- — hO *t-$ , (17) 

/$$rfR ;$^ T rfR j^t^t^r 
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where X is a Hermitian operator. This expression can be obtained by Taylor expansion 
in (\& T - $). 

Application of equation ( TlTI) to the nodal term of equation ( fl6l) gives 



IW T d^T ,o 

1 7s ^t^t dA 

2 ^$$rfR 



t y d^T 

^ T ^ T dA 



$\[> T dR 



T T W T d^ T JO 
S ^ T ^T dA 



+ 

2 I ^ T ^ T (iR 

+ 0[(^ T - $) 2 ] . (18) 

It is straightforward to show that the variational term in equation (1181) is zero because 
is continuous and different iable across the nodal surface [35]. The mixed-distribution 
term on the right hand side of equation (1181) is, moreover, exactly equal to twice the 
second term on the right hand side of equation (fTTT) [35] , and therefore 



W T d ^T ,„ 

1 7s ^ T ^T dA 



-E D ) d^ T 

dT 



dR 



$$dR y^T^rfR 

+ 0[(* T -$) 2 ]. (19) 

The mixed-distribution term on the right hand side of equation ( 1191) may be evaluated 
straightforwardly if H is taken to act to the right. Next we apply the extrapolation 
procedure again to eliminate the mixed distribution term in equation ( fl9l) . which leads 
to 



_ , V*t d* T JO 
1 7s grj^T dA 
f $$dR 



(H — Ed) d^ 



dA 



dR 



+ 



$$dR 



dA 



^T^dR 



0[(* T - $) 2 ] . (20) 



Using equations (]14p . (I16p . and (j20p with H acting to the left in the variational term, 
we obtain our final expression for the energy derivative [31] 



dE D . 




"1 dH 
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$ dA 
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(21) 
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Equation ( 121 j) consists of integrals over the pure and variational distributions. This 
expression is expected to give much more accurate forces than equation pip with 
the approximation ( Tl2l) because (i) the Pulay terms with the pure distribution are 
smaller than those obtained with the mixed distribution, and (it) the error is of second 
order rather than first order. Equation (121j) is expected to show reasonable statistical 
behaviour because it satisfies a zero- variance principle. If \I/t and d^x/dA are exact, 
the pure distribution term in equation (121]) gives the exact result even for a single 
configuration, and the variational term has the same property and, in addition, gives 
zero if \I/t is exact. 

If the nodal surface of \Px is independent of the variable parameter in the 
Hamiltonian then only the Hellmann-Feynman term survives in equation (l2~Tj) . The 
fixed-node approximation is equivalent to placing infinite potential barriers everywhere 
on the nodal surface of \I/x- The nodal surface term can be viewed as arising from 
the change in the potential barriers as the variable parameter changes. Alternatively, 
one can include the infinite potential barriers in the Hamiltonian, formulating a fixed- 
node Hamiltonian Hp^ [35]. The Hellmann-Feynman theorem holds for Hfn, but iJpN 
depends on the nodal surface of \I/x and the derivative dHp-^/dX generates the same 
Pulay terms as in the approach described above. 



5. Total versus partial derivatives of \l/x 

We assume that the trial wave function \l/x (and hence the nodal surface) depends on 
N c variational parameters q, with % = 1, ...,N C . We further assume that \l/x has both 
explicit and implicit dependence on the parameter A in the Hamiltonian we wish to 
vary. For example, the atomic centres of a localised basis set or the Jastrow factor may 
have explicit dependence on the positions of the atoms, and they may contain variable 
parameters whose optimal values depend implicitly on the atomic positions. The total 
force calculated with the variational, mixed or pure distributions can be written as 



dE dE Sk dE dc- 



dA 



A=0 



a=o ~ dd dA 



(22) 

A=0 



E 

8=1 

where E = Ey for VMC and E = E^> for the mixed and pure DMC distributions. 
The two terms on the right hand side of equation ( 122]) correspond to the explicit and 
implicit dependences on A. The choice of these dependences is not unique and may 
depend on how the parameters are considered in a QMC simulation. It is possible to 
construct \I/x without variational parameters if the values of all of the parameters q 
are kept fixed. Also, \I/x can be constructed so that it has no explicit dependence on 
the nuclear coordinates when, for example, the centres of atomic-centred basis sets are 
made variational parameters [26]. For QMC calculations it is, however, convenient to 
choose forms of \l/x where both explicit and implicit dependences on A are present. 

Calculating the force using equation (1221) involves knowledge of how the wave 
function and hence the q change with A, i.e., dcij/dA. In standard quantum chemistry 
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methods, these derivatives are typically obtained by second-order perturbation 
theory [36]. Unfortunately such an approach is not straightforward in VMC and DMC. 

We follow a different route and approximate all total derivatives by partial 
derivatives or, equivalently, the term involving the sum in equation (1221) is neglected. 
This approximation introduces an error in the force of first order in (\I/t — We expect 
this approximation to be rather accurate in both VMC and DMC for the following 
reasons. In VMC each term in the sum vanishes if E is energy-minimized with respect 
to the corresponding q, i.e., dE/dci = [26]. In practice this condition will not normally 
be satisfied exactly because i) the q are obtained by stochastic methods and are therefore 
subject to statistical noise, ii) the q are sometimes obtained by minimizing the variance 
of the local energy rather than the energy itself, in) the orbital parameters are often 
obtained from density-functional theory calculations and are therefore energy minimised 
with respect to the wrong Hamiltonian. Similar considerations apply within DMC, and 
we expect that E is approximately minimized with respect to the q so that the dE/dci 
are very small. The terms in the sum in equation (1221) are therefore expected to be very 
small and can be neglected. 



6. Practical test of DMC forces 



Before presenting some results, it is worth discussing the accuracy we require of QMC 
forces for them to be useful. It is helpful to think about the error in the force on an atom 
in a molecule or solid in terms of the associated error in the bond length. The errors in 
equilibrium bond lengths obtained from well-converged density functional theory (DFT) 
calculations for sp-bonded atoms are (odft — o exp ) ~ 0.01 A, where "exp" denotes the 
experimental value. To be useful, the errors in DMC bond lengths measured as the 
difference between the values obtained from the forces and energies (admc 3 — a DMcf eS ) 
should be substantially smaller than ~0.01 A. 

We illustrate the results obtained for the mixed and pure DMC force estimates 
with calculations on the SiH molecule. It should be noted that we have calculated the 
forces on both the Si and H atoms, which should of course be of equal magnitude and 
opposite sign. This is in fact a highly non-trivial test of the computational method as 
the condition of zero net force on the molecule is not automatically satisfied in these 
calculations. It is easier to obtain accurate trial wave functions for H atoms than for 
heavier atoms and therefore the forces on the H atoms tend to be more accurate. A 
number of other DMC studies have calculated forces on small molecules containing H 
atoms in which the bond lengths were deduced entirely from the forces on the H atoms. 
We do not regard this procedure as a proper test of the methodology. 

In our calculations we used Dirac-Fock based non-local pseudopotentials to 
represent the H + and Si 4+ ions [371 1381 139] . The Hellmann-Feynman contributions 
to the forces from the non-local pseudopotentials were calculated using the expressions 
given in the Appendix of Ref. [29]. The pseudopotential localisation procedure [10] 
introduces additional terms in the force which contain d\&T / dA [41] . These terms arise 
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Figure 1. The mixed estimates of the Hcllmann-Feynman force (F^f^) and the total 
force (F^J evaluated for the H and Si atoms. The derivative of the Morse potential 
fitted to the DMC energies (d-En/dA) is shown as reference. 

from the derivative of the Hamiltonian and therefore we have included them in the HFT 
expression for the forces in the following discussion. 

The trial wave functions were constructed as follows. We calculated the molecular 
orbitals within Hartree-Fock theory using the GAMESS code [12] with a fairly small 
basis set consisting of 5s and 2p Gaussian basis functions for each atom. We used 
a Jastrow factor which included electron-electron and electron- ion terms jl3]. The 
initial values of the Jastrow parameters were obtained by minimising the variance of the 
energy [H] and they were further refined by minimising the VMC energy |l5j 06j H7] . 
We could have generated considerably more accurate trial wave functions by using a 
larger Gaussian basis set, a more complicated form of Jastrow factor, including several 
determinants, or backflow transformations jH]. We chose not to do this because we 
wanted to work with trial wave functions of a similar quality to those we could readily 
generate for heavier atoms and much larger systems. When evaluating d^x/dA we 
included only the explicit dependence on the ionic positions which occurs in the Jastrow 
factor and Gaussian basis functions. Expectation values with the pure distribution were 
evaluated using the future- walking method [23] . 

Results using the mixed distribution expression of equation ( TTTT) with the 
approximation of ( JT2I) are plotted in figure [TJ As reference, we have also plotted the 
derivative of the Morse potential previously obtained from a fit to the DMC energies. 
The minimum in the DMC energy occurs at a bond length of 1.5242(6) A. The forces 
on the H and Si atoms should be equal and opposite, and figure [T] shows that the 
total forces F^ on the H and Si atoms sum almost to zero over the range of bond 
lengths shown, with an error corresponding to an error in the bond length of equal or 
less than 0.001 A. However, the total force gives an equilibrium bond length which is 
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Figure 2. The same as figure HJ but for the pure forces evaluated using equation (|21|) . 

about 0.010 A too short. The Hellmann-Feynman force F^f^ on the H atom is zero at 
a bond length which is 0.003(1) A too long, while that on the Si atom is zero at a bond 
length which is 0.020(1) A too short. Including the Pulay terms improves the agreement 
between the forces calculated on the H and Si atoms, but overall the bond lengths from 
the mixed-estimator forces are no better than those which could be obtained in DFT 
calculations. 

Results using equation (l2"Tj) are plotted in figure |2J We refer to this as the expression 
for the pure forces, although it also contains a term evaluated with the variational 
distribution. The total pure forces F^ rc on the H and Si atoms sum almost to zero 
over the range of bond lengths shown, with an error corresponding to an error in the 
bond length of less than 0.003 A, which is about three times the statistical error. The 
equilibrium bond lengths from the zero force condition are also within statistical errors 
of the equilibrium bond length obtained from the minimum in the DMC energy. The 
nodal terms in equation f[2"Tj) are by no means negligible. The pure Hellmann-Feynman 
forces -F^rc" on the H and Si atoms do not sum to zero, and the force on the H atom 
is zero at a bond length which is 0.008(1) A too long, while the force on the Si atom 
is zero at a bond length which is 0.010(2) A too long. The total pure forces are much 
more accurate than the total mixed forces, which suggests that the approximate forms 
of equation (l2~Tj) and d^-p/dA give rather accurate results in our simulations. The bond 
lengths obtained from the total forces and from the HFT term, and using forces on the 
Si and H atoms, within mixed and pure DMC, are summarised in Table [TJ It is clear 
that the total pure forces give equilibrium bond lengths differing from that obtained 
from the DMC energies by much less than the typical error in DFT bond lengths of 
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7. Statistical properties of the force estimators 

A statistical estimate of a quantity should be accompanied by a confidence interval which 
indicates the reliability of the estimate. Consider a probability distribution function 
(PDF) P(x) which has a mean value \x and variance a 2 . The mean value of a random 
variable x drawn from P(x) may be estimated using the sample mean, 

1 N 

X N = J-jY, X ii ( 23 ) 
iV t=l 

and the variance estimated using the sample variance, 

1 N 

a * = jv^iS (Xi "^ )2 ' (24) 

with a 



We then say that xn lies within the interval xn — o~n/ v N , x~n + a^/v N 
confidence of 68 %. 

This statement relies upon the probability distribution satisfying certain conditions 
which are the content of the Law of Large Numbers (LLN) and the Central Limit 
Theorem (CLT). The LLN concerns the statistical convergence of the sum in equation 
( 123]) to the true mean value with increasing N. Similarly, the CLT concerns the 
statistical convergence with increasing N of the probability distribution from which 
xn is drawn to a Normal distribution with mean \i = fxP(x)dx and variance 
a 2 /N = l/N J(x — /i) 2 P(x) dx. For the LLN to be valid it is sufficient that the first 
moment of P(x) exists, which requires that P(x) decays faster than \x — xq\~ 2 , where 
Xq is a constant. Similarly, the CLT is valid if the second moment of P(x) exists, which 
requires that P(x) decays faster than \x — xo| -3 - 

Equations (JSJ) and (JS} show that the VMC and DMC energies are evaluated as 
averages of local energies E^{Hi). The probability distribution of the local energies 

| By this we mean that the convergence is "almost sure" in the sense that the underlying probability 
density function approaches a delta function in the limit of a large sample size. 



Table 1. Equilibrium bond lengths a in A for the SiH molecule obtained from the 
total forces and from the HFT term and using forces on the Si and H atoms, within 
mixed and pure DMC. The equilibrium bond length obtained from the minimum in 
the DMC energy curve is 1.5242(6) A. 





mixed DMC 


pure DMC 


a tot (Si) 


1.5140(3) 


1.5259(8) 


a HFT (Si) 


1.5044(8) 


1.5343(16) 


a tot (H) 


1.5150(3) 


1.5227(7) 


a HFT (H) 


1.5272(9) 


1.5320(12) 
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-P(-El) is generally non-Gaussian and has "fat tails". The asymptotic behaviour of 
-P(-El) is determined by the singularities in E^. Trail [49] considered singularities in E-^ 
arising from a Coulomb energy divergence at particle coalescences, divergences in the 
local kinetic energy at the nodal surface of ^t, and in finite systems where \I/t has the 
wrong asymptotic behaviour. The Coulomb energy divergence at particle coalescences 
can be removed by imposing the Kato cusp conditions on \&t> as is normally done 
in QMC calculations, but it is not clear how to remove the divergences at the nodal 
surface. Trail [19] showed that the singularities in E^ lead to tails in P(Ei) which 
decay as \E L — E \~ 4 , where E is a constant. Consequently the energy estimate is 
drawn from a Normal distribution whose variance we may estimate (in the limit of large 
N). However, the estimated variance of the local energy is drawn from a distribution 
of infinite variance that is not Normal, hence it is not entirely clear how to obtain a 
meaningful confidence interval for this estimate. 

A similar analysis can be performed for the forces. Each term in the force can be 
written as a sum of local contributions, 
1 N 

^^E^(R*). (25) 

i=l 

For the Hellmann-Feynman term we distinguish between three cases: an all-electron 
atom, a smooth local pseudopotential, and a smooth non-local pseudopotential. The 
probability distribution of the Hellmann-Feynman term for an all-electron atom has 
slowly decaying tails of the form \F^ — F |~ 5 / 2 , where F is a constant. In this case the 
LLN applies but the CLT is not satisfied, and the variance of the Hellmann-Feynman 
force is infinite. Although a sample variance can be calculated it does not statistically 
converge to a constant value with increasing N, and it is not related to the width of 
the distribution from which the estimated force is drawn. The infinite variance in the 
Hellmann-Feynman force for an all-electron atom has long been recognised and methods 
for dealing with it have been devised. Assaraf and Caffarel [251 12Z] suggested adding 
a term to the HFT force which has zero mean but may greatly reduce the sample 
variance. (This amounts to adding the term containing d^x/dA in the first term on 
the right hand side of equations (JTUJ) and f[2"T|) .) A similar approach was developed 
by Per et al [50]. Chiesa et al [51] developed a method to filter out the component 
of the electronic charge density responsible for the infinite variance. Using smooth 
local or non-local pseudopotentials without singularities at the nucleus also eliminates 
the infinite-variance problem arising from the electron-nucleus interaction. The terms 
in the force estimators which contain d\I/x/dA also give probability distributions with 
I-^l ~~ -^o|~ 5 ^ 2 tails, and hence the CLT is not applicable. The origin of the problem is 
that although \&x is zero on the nodal surface, d\I/x/dA is normally finite. Trail [52] and 
Attaccalite and Sorella [30] have devised rigorous methods within VMC to eliminate the 
infinite variance by altering the probability distribution near the nodes. 

Figure |3] shows probability distributions for VMC local forces on the Ge atom in a 
GeH molecule. The upper graph shows that the probability distribution for the sum of 
the local Hellmann-Feynman force (the term in equation ( flOl) containing dH/dX) and 
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Figure 3. The probability distribution of the VMC local Hellmann-Feynman forces 
(red) and the local Hellmann-Feynman forces with the zero variance force (blue). 
The forces on the Ge atom of a GeH molecule are shown calculated using non-local 
pseudopotentials and a Slater- Jastrow wave function. The distributions were generated 
from 10 7 samples. The upper and lower graphs show the same data on different scales. 

the local zero-variance force (the term in equation ( fTOj) containing H) is much more 
strongly peaked in the region of the mean value than the local Hellmann-Feynman force 
itself. This illustrates the operation of the zero-variance term. The lower graphs show 
the same data plotted on different scales. The inset graph shows a very wide scale of 
local forces, and the local zero-variance force actually has more outlying contributions 
than the local Hellmann-Feynman force. This illustrates the result that the probability 
distribution of the local Hellmann-Feynman force has \Fl — Fq\~ 4 tails (for a smooth 
pseudopotential) and the CLT holds, while the probability distribution of the local 
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zero- variance force has — F \~ 5 / 2 tails so that the CLT is not satisfied. This seems 
somewhat counterintuitive; without the zero-variance term the variance is finite but 
does not approach zero when (^t, d^x/dA) are exact whereas, if we include the zero- 
variance term, the variance is zero for the exact case, but infinite for the approximate 
case. Analogous results apply for the DMC forces obtained with the mixed and pure 
distributions. 

Currently no rigorous procedure has been demonstrated for removing the infinite 
variance which arises in DMC calculations of the Pulay forces. Although this state of 
affairs is unsatisfactory, the probability distribution of the local forces is strongly peaked 
in the region of the true value, which suggests that reasonable procedures for estimating 
a confidence interval can be devised. 

8. Conclusions 

We have formulated DMC forces evaluated with the pure distribution. The force can 
be written as the sum of a Hellmann-Feynman force and a term which can be cast as 
an integral over the nodal surface. The nodal term can be recast as an average over the 
pure distribution of a quantity which depends only on \I/x and d\I/x/dA. It is obvious 
that such an exact expression must exist, as \1/t and d\I/x/dA fix both the nodal surface 
and its first order variation, which in turn fix the DMC energy and its first derivative. 
It is however, very awkward to evaluate an integral over the nodal surface, and therefore 
we have adopted an approximate approach to evaluating it which has an error of order 
(\1/t — $) 2 . Our final expression ( 121]) consists of integrals over the pure and variational 
distributions, both of which exhibit a zero variance property. 

The forces evaluated using equation (j2Tl) yield bond lengths for the SiH molecule 
with an error (with respect to DMC energy calculations) of much less than 0.01 A, 
which is acceptable for most purposes. Of course, before claiming that the problem 
of calculating DMC forces is "solved", accurate results for heavier atoms and larger 
systems are required. To this we should add the issues of efficient generation of the pure 
probability distribution for large systems, overcoming the "infinite variance" described 
in Section[71 obtaining accurate forms for \l/x and d^x/clA, and (hopefully) removing the 
second approximation in equation fl21]) . Notwithstanding these caveats, we believe that 
substantial progress has been made in calculating DMC forces, and that their evaluation 
will become routine in the not-too-distant future. 
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